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SUMMARY 

Highly convective scalar transport involving near-discontinuities and strong streamline 
curvature was addressed in a paper by Smith and Hutton in 1982, comparing several differ- 
ent convection schemes applied to a specially devised test problem. First-order methods 
showed significant artificial diffusion, whereas higher-order methods gave less smearing but 
had a tendency to overshoot and oscillate. Perhaps because unphysical oscillations are more 
obvious than unphysical smearing, the intervening period has seen a rise in popularity of 
low-order artificially diffusive schemes, especially in the numerical heat-transfer industry. 
The present report describes an alternative strategy of using non-artificially diffusive higher- 
order methods, while maintaining strictly monotonic transitions through the use of simple 
flux-limiter constraints. Limited third-order upwinding is usually found to be the most cost- 
effective basic convection scheme. Tighter resolution of discontinuities can be obtained at 
little additional cost by using automatic adaptive stencil expansion to higher order in local 
regions, as needed. 


1. INTRODUCTION 

In a well-known paper published in 1982, Smith and Hutton presented results of several 
authors’ attempts to numerically solve a specially devised test problem involving streamline 
curvature typical of recirculating flows and steep variations in the transported scalar 1 . Most 
schemes were able to handle the diffusion-dominated low-Peclet-number regime adequately; 
but in the important high-convection regime, Smith and Hutton concluded that convection 
modelling “remains the art of compromise between diffusive and oscillatory errors.” In the 
intervening period, it seems that artifically diffusive low-order (blended first/ second-order) 
convection schemes have become more popular than higher-order potentially oscillatory 
methods. Perhaps this is because the overshoot (or undershoot) problems associated with 
higher-order methods can lead to obviously unphysical results such as locally negative 
densities or absolute temperatures, for example. But the low-order methods’ results are also 
usually highly unphysical — although perhaps not always obviously so. 

Convection-diffusion schemes that revert to first-order upwinding for convection (while 
physical diffusion is ignored) under high-convection conditions achieve (plausible looking) 
nonoscillatory results by replacing the high-convection physical problem, with an artifically 


♦Work funded by Space Act Agreement C-99066-G. 



diffusive (and anisotropic) /ow-convection numerical problem. Blended flrst/second-order 
schemes of this type, such as Spalding’s “Hybrid” method 2 and Patankar’s “Power-Law 
Difference Scheme” (PLDS) 3 seem to have gained in popularity in recent years, especially 
in the numerical convective heat-transfer industry 4 . This is puzzling, because there have 
been several studies showing the artificially diffusive nature of such schemes 5 ' 10 over the past 
decade or more. Even more puzzling is the fact that many authors use Hybrid, PLDS or 
similar exponential-based schemes" in combination with sophisticated (and expensive) 
multiple-equation turbulence models — apparently not realizing that the turbulence model 
is being used throughout most of the flow-field merely as an expensive diagnostic to switch- 
off the physical (turbulent plus laminar) diffusion terms in the governing equations, replacing 
them with (anisotropic) artificial numerical diffusion. It is perhaps not surprising that the 
results of such computations are typically very insensitive to the turbulence model being 
used. Most puzzling of all is that some researchers continue to use first-order-based 
convection schemes to actually develop and “tune” new turbulence models. 

The usual “justification” for this approach seems to be based on grid-refinement 
studies, i.e., the grid is refined to a point where the results do not seem to be changing very 
much. But, for first-order methods, the approach to true grid-independence is a notoriously 
slow process. One cannot claim reasonable accuracy (or proper use of the turbulence model) 
until the grid is refined to a point where the component grid Peclet (or Reynolds) numbers 
are everywhere 0(1) or less — in which case Hybrid and PLDS are operating as second- 
order central differencing. The massive grid refinement that this would call for under high- 
convection conditions is clearly impracticable. 

There is a strong need for a conceptually simple (and computationally inexpensive) con- 
vection scheme giving highly accurate non-artificially diffusive and non-oscillatory results 
on practical grids under high-convection conditions; i.e., for grid Peclet (or Reynolds) 
numbers arbitrarily large. The same scheme should, of course, be able to handle the low- 
Peclet-number diffusion-dominated regime, as well. As will be shown in this report, these 
apparently conflicting requirements are not incompatible. Third-order upwinding is the 
lowest-order convection scheme for which the leading truncation error is dissipative (involv- 
ing even-order spatial derivatives) but not “diffusive” (i.e., second-order derivatives) — by 
definition, leading truncation error in this case involves fourth-order spatial derivatives. As 
is well known, however, third-order upwinding in its basic form can give rise to unphysical 
overshoots or undershoots near regions involving rapid changes in the transported variable. 

But it is a relatively simple matter to incorporate universal limiter constraints (appli- 
cable to any order of accuracy) giving tight monotonic resolution of near-discontinuities 
without corrupting the accuracy of the underlying scheme. This universal limiter for tight 
resolution and accuracy implemented via a simple high-accuracy resolution program consti- 
tutes the ULTRA-SHARP strategy for high-convection modelling. The recommended 
method uses limited third-order upwinding (ULTRA-QUICK) as the basic convection- 
diffusion scheme; then, in local regions requiring even higher-order resolution, the algorithm 
automatically branches to a limited higher-order scheme (ULTRA-5th or ULTRA-7th 
upwind, for example) using adaptive stencil expansion, locally, controlled by a simple non- 
smoothness monitor. In terms of achieving a desired accuracy (compared with a known 
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exact solution, for example), this strategy is optimal in terms of requiring the lowest 
computer cost as the grid is refined. In other words, although low-order methods are less 
expensive per grid point, they require an exorbitantly fine grid to achieve a prescribed 
accuracy; by contrast, the higher cost (per grid point) of very high-order methods used 
globally is not completely offset by the lower cost of a concomitantly coarser grid. 

In the following sections the Smith-Hutton test problem is briefly described. Then, for 
reference, results are shown for a number of well-known convection-diffusion schemes; in 
particular: PLDS (representative of exponential -based schemes), second-order upwinding, 
and third-order upwinding (QUICK), using a conservative control-volume time-marching 
formulation. Fifth- and seventh-order upwinding are also briefly discussed; as is typical of 
unlimited higher-order schemes, tighter resolution in this case is offset by stronger 
oscillations. The concept of the universal limiter, based on normalized variables, is then 
briefly reviewed. Results are shown for limited third-order (ULTRA-QUICK) and an 
ULTRA-3rd/5th/7th-order scheme using local adaptive stencil expansion. Finally, a cost- 
effectiveness study shows the optimality of the third-order-based ULTRA-SHARP schemes. 

2. TEST PROBLEM 

The two-dimensional test problem devised by Smith and Hutton is concerned with steady- 
state convection and diffusion of a scalar field such as temperature, T, for example, in a 
prescribed velocity field, v(x, y), with a known constant diffusivity, D. The nondimen- 
sional governing equation is 


v-V T = — V 2 T 
Pe 


( 1 ) 


introducing the (macroscopic) Peclet number 

Pe = — ~ = const ( 2 ) 

D 


using appropriate reference velocity and length scales. The flow domain considered is a 
rectangle: - 1 <jc < 1 , 0<y<l. And the velocity field is given by 

u = 2y (1 - x 2 ) (3) 


and 


v = - 2x (1 - y 2 ) 


(4) 
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corresponding to a streamfunction 


* = - (1 - x 2 ) (1 - y 2 ) 


(5) 


Figure 1 shows the streamline pattern for this flow-field. 

The inlet temperature profile is specified as 

TJx ) = 1 + tanh[a(l + 2x )] (6) 


for y = 0 and -l<;t<0. This is also shown in Figure 1. For x = -1, the left-hand 
boundary condition becomes 


T„ * TJ-1) = 1 - tanh « (7) 

This is used as the boundary condition along the boundary streamline \p b = 0; i.e., at 
x — ± 1 (for 0<y<l) and at y = 1 (for -1<*<1). For a greater than about 3, this 
means that the boundary temperature is essentially zero, whereas the top of the inlet profile 
is very close to 2 as x -*> 0. Smith and Hutton proposed a = 10 as representative of a 
relatively sharp transition. In the present paper, two other values of a are used: a = 100 
(representative of a very sharp transition) and a = 5 (representing a relatively smooth 
transition). Note that no physical boundary conditions are specified at the outlet boundary, 
y = 0 (0<jc< 1). Numerical boundary conditions equivalent to 

= 0 for 0 < jc < 1 (8) 

y-Q 


d T 
Ty 


are described below. 

In the case of purely convective flow, Pe -*• <», the exact solution is easily obtained, 
since T = const along streamlines; i.e., T = T For example, at the inlet, 

x = - JTTJ , - 1 < jc < 0 (9) 

so, throughout the flow domain when there is no diffusion, 

T ($) = 1 + tanh[a(l - 2 \fTTJ )] = T(x, y) (10) 
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using Equation (5). In particular, this gives an outlet profile as the mirror-image of the inlet 
profile 


F out (x) = 1 + tanh[a(l - 2x)] 


(ID 


for y = 0 and 0 < x £ 1. Figure 2 shows a three-dimensional portrayal of Equation (10) 
on a 40x20 uniform mesh (41 X 21 grid -points, with Ax = Ay). In Figure 2(a), a = 100; 
whereas in 2(b), a = 5. 

Figure 3 shows portions of a typical staggered mesh used in the present analysis. Note 
that F-nodes are placed at boundaries. Boundary nodes shown as dots within squares 
correspond to specified boundary conditions; solid dots represent interior computed F-nodes; 
exterior pseudo- F-nodes (triangles) are also shown for use with higher order methods. 
Hollow circles represent ^-nodes; these occur at the corners of temperature control-volume 
cells. This is shown in more detail in Figure 4. This is a convenient arrangement, since 
average cell-face convecting velocities are then available by simple subtraction of stream- 
function values; e.g., referring to Figure 4, 


u t 



( 12 ) 


so that the average left-face Courant number is 


CXL(i, j) 


WtL ^ BL ) ^ 

Ax Ay 


(13) 


Computation of interior node values of F follows a simple time-marching procedure. 
One first computes the left and bottom face fluxes 


and 


FLUXL(i, j) 


CXL • T t 


CXL| Ax 
PXL 


d T 
d x t 


(14) 


FLUXB(i, j) 


CYB • T b 


CYB| Ay 
PYB 



(15) 
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introducing local component cell face Peclet numbers 

PXL = |u, | Pe Ax 
and 

PYB = |vj Pe Ay 

The new values of T are then updated by a simple assignment statement, 

Set: T(i,j) = T(i, f) + FLUXL(i, j) - FLUXL(i+l, j) 

+ FLUXB(i, j) - FLUXB (i, j+ 1) 

where flux conservation has been observed at each face. This is repeated until a converged 
steady state has been achieved. 

The treatment of the outflow numerical boundary condition is shown in Figure 5. 
Assume that T (y) follows a parabola near the outlet for y ^ 0, with three conditions 


T(2Ay) = T 2 

(19) 


T(Ay) = T x 

(20) 

and 

' d t 
d y 

K J 

= 0 

0 

(21) 

Then, it is not hard to show that the corresponding value of T at the boundary is 


T 

4 

3 

T. - - T, 
1 3 2 

(22) 

The corresponding pseudonode values of T. 
upwinding) are taken to be simply 

i and T -2 (used for fifth- and seventh-order 

T’-i 

= ! 

r* 

to 

II 

(23) 


( 16 ) 

(17) 

(18) 


as shown in the figure. Pseudonode values along the other boundaries are set equal to their 
adjacent boundary values in a normal direction. 
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3. EXPONENTIAL-BASED SCHEMES 


Exponential-based convection-diffusion schemes were first introduced into computational 
fluid dynamics by Allen and Southwell 12 , and have been “rediscovered” in various 
equivalent or approximate formulations by several people in the past thirty-five years. 
Spalding’s Hybrid scheme 2 , Patankar’s PLDS 3 , and the algebraic approximation of Raithby 
and Schneider 11 can be interpreted as various levels of approximation to the exponential 
differencing scheme (EDS). It is fairly easy to show 13 that EDS is equivalent to using 
second-order central differencing for both convective and diffusive fluxes while replacing 
the actual grid Peclet (or Reynolds) number, P A , with an effective value, P A \ that is itself 
a function of the physical P A . The functional relationship is 13 

P A * = 2 tanh(P A /2) ( 24 ) 

Spalding’s Hybrid scheme can be interpreted as a very rough approximation to this, given 
by 


= P. 


for 0 < P A < 2 


for P A > 2 


(25) 


Patankar’s power-law difference scheme represents a much more accurate approximation of 
the hyperbolic-tangent function 


2P a 


P A + 2(1 - 0.1 P A ) 5 


for 0 < P A < 10 


p; - 2 


for P A > 10 


(26) 


The algebraic formulation of Raithby and Schneider can be interpreted as 

P* = P [ + (1 + Q 005 1 1 

2(5 + P A 2 ) P A (1 + 0.05 P 2 ) 


(27) 


Note that here, too, P A -* 2 for large values of P A . In fact, for EDS itself, Equation (24), 
P A * 2 for P A > 6 (since tanh 3 = 0.995... ). All three approximations are shown in 
relation to the hyperbolic-tangent curve in Figure 6. 
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For exponential-based schemes, the left-face flux, for example, is given by 

FLuxL(i, i) = £xl (j t j * r Mi/ ) - -^iL (T u - r M ,) (28) 

where PXL' is the effective local ^-component grid Peclet number at the left face. Note 
that for PXL* = 2, the flux becomes 

FLUXL(/, f) = CXL*r_ i ; for CXL > 0 (29) 

or 

FLUXL(i, j) = CXL -7^ for CXL < 0 (30) 

This, of course, corresponds to first-order upwinding for convection, with physical diffusion 
(computed but) ignored. This is what occurs in the Hybrid scheme for P A > 2, and in the 
other schemes (including EDS) for P A greater than about 6. 

Figure 7 shows 40x20 results for Pe = oo ; in this case the scheme is operating 
everywhere as first-order upwinding. By comparison with Figure 2, these are seen to be 
very artificially diffusive results. This is typical of exponential-based schemes. Figure 8 
gives inlet and outlet profiles using a = 100 for Pe = oo , 500, and 10, showing com- 
puted solutions on 20 x 10, 40x20, and 80x40 grids in each case, using PLDS. The refer- 

ence finite-Peclet-number results have been obtained using the ULTRA-3/5/7 upwind scheme 
(described later) on a very fine (160x80) grid. For the larger-Pe cases, the gross artificial 
diffusion of the exponential-based scheme is clear. In the case of Pe = 10, local grid Peclet 
numbers are small and the scheme is equivalent to second-order central differencing. The 
error reported in the figure captions is computed using 

* - ^ El r ~, - T ,J < 31 > 


where the summation is over all interior grid-points plus the outlet boundary, and N is the 
total number of grid-points involved, excluding pseudonodes (i.e., N = 21 x 11, 41 x21, or 
81x41). 
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4. SECOND-ORDER UP WINDING 

When second-order upwinding is used for convection, it is conventional to use second-order 
central differencing for diffusion. In this case, the left-face flux, for example, is given by 

FLUXL (i,j) = CXL-r, - (T u - T,_, ,) (32) 

where the left-face convected value is 

T, = 1 T mj - 1 r 2j for CXL > 0 (33) 

or 

T, = | T., - -i T M j for CXL < 0 (34) 

Similar formulas are easily obtained for the bottom-face flux. Equations (33) and (34) can 
be combined into a single form valid for positive and negative convecting velocities by 
writing 

T, = - (T t , + T tl ,) - I CURVNL (35) 

r 2 J 1-1 2 

defining the “normal curvature” at the left face as 

CURVNL = CRVAVL - SGN(C ^H THIRDL (36) 

2 

where (suppressing the y'-index, for convenience) the average (symmetric) second-difference 
across the left face is 

CRVAVL = I (7;. +1 - + r ( . a ) (37) 
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and the third-difference across the face is 


THIRDL = T itl -3 T ( + 37;.., - T., 2 


(38) 


The infinite-Pe results for second-order upwinding are shown in Figure 9 — as usual, 
for a = 100 and 5. In this case, the smooth-inlet-transition results are quite good, with 
only a little numerical spreading and a very slight overshoot near the outlet; but note the 
significant overshoots and undershoots in the sharp-transition case. Figure 10 shows grid- 
dependence results for Pe = <» , 500, and 10, for a = 100. 

5. THIRD-ORDER UPWINDING (QUICK) 

The QUICK scheme (quadratic upstream interpolation for convective kinematics) is the 
canonical third-order-upwind scheme for steady-state flow 14 . In this case, the left-face flux, 
for example, has the same form as Equation (32); however, for consistency, the convected 
face value includes both normal and transverse curvature effects 14 and the normal-curvature 
coefficient is much smaller 

T t = 1 (T t J + T MtJ ) - i CURVNL + .-L CURVTL (39) 

where CURVNL is given by Equation (36) and the upwind-weighted transverse-curvature 
term is 

CURVTL = r M>y+I - 2T m j + 7;.., .., for CXL > 0 (40) 

or 

CURVTL = r (t/4l - 2T. j + T. for CXL < 0 (41) 

Note that, consistent with bi-quadratic interpolation in the vicinity of the left face, the 
diffusive flux is identical to that obtained by central -differencing 14 . 

The QUICK results for Pe = °° are shown in Figure 11. The smooth-transition case 
is very well modelled; but, as with second-order upwinding, overshoots and undershoots 
occur in the sharp-transition simulation, although the computed transition itself is noticeably 
sharper in this case. Grid-dependence results are seen in Figure 12. 
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6. FIFTH- AND SEVENTH-ORDER UPWINDING 

The fifth-order upwind algorithm used in this paper again takes the form of Equation (32), 
but in this case 

T ( = i (T l } + T^ j) - i CRVAVL + FORTHL + A. CURVTL (42) 

where (again suppressing j's for convenience) the upwind- weighted fourth difference is 
FORTHL = T ; +1 - 47;. + 6T ' M - 47;._ 2 + T._ 3 for CXL > 0 (43) 

or 

FORTHL = T M - 4 T M + 6T t - 47].., + T t _ 2 for CXL < 0 (44) 


Three points should be mentioned: 

(i) Higher-order terms are not used in the diffusive flux. This is appropriate 
because, when diffusion is large (small Pe), modelled profiles are smooth and 
the second-order form is entirely adequate; whereas, under high-convection 
conditions, the form of the small diffusion terms is not very important. 

(ii) Higher-order transverse terms are not used in the convective flux. Although the 
third-order transverse curvature term is significant, numerical experimentation 
has shown that higher-order transverse terms have an almost negligible effect on 
results; but inclusion would add significantly to the cost of the calculation. 

(iii) The coefficient of the normal-curvature term (1/6, rather than the theoretical 
value of 1/8) has been found to give slightly more accurate results in cases of 
scalar convection and diffusion, where exact solutions are available 13 . This was 
not found to be the case with third-order upwinding — where 1/8 seems to be 
optimal in all cases tested. 

Figure 13 gives the fifth-order results for Pe = oo . As perhaps expected, the large-a 
transition is sharper (than third-order) but generates significantly more overshoots, under- 
shoots, and secondary ripples. The smooth transition is graphically indistinguishable from 
the exact result. Grid-dependence studies are again predictable and need not be shown here. 
Higher (for example, seventh) order upwinding merely accentuates the trends seen with fifth- 
order. 
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The seventh-order formula used in this study takes the form 


T t = - (T i } + r.., y ) - 1 CRVAVL + FTHAVL 


1 SIXTHL + — CURVTL 


100 


6 

24 


128 


( 45 ) 


where (suppressing j’s, as usual) the average (symmetric) fourth-difference across the left- 
face is 


FTHAVL - i (T hl - 3T„, * IT, * 2T M - 3T,_ 2 + T M ) 


(46) 


and SIXTHL is the upwind-weighted sixth-difference 


SIXTHL 


T 

x i* 2 


- 6 T„, * 15 T, - 207;.., * 157- m - 67;.., * 


T 


i-4 


(47) 


for CXL > 0; all indexes in SIXTHL are increased by 1 for CXL < 0. 

7. UNIVERSAL LIMITER 

The universal limiter is most easily described in terms of normalized variables. Let T f 
represent the value of the convected scalar at any control-volume face; call the adjacent 
downstream node value T D , the adjacent upstream node value T c , and the next upstream 
node value (in a direction normal to the face) T v . Figure 14 sketches the definition of these 
terms; as seen, node C lies between nodes U and D . Note, however, that the nodes in- 
volved are dependent on the sign of the normal convecting velocity component, u n , at the 
CV face. Now define, anywhere in the vicinity of the face, a normalized variable 


T(x, y ) 


T(x, y ) - T v 
T D ~ T v 


(48) 


In particular, 


T f 



(49) 
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and 


T - T 

= 'c *_u (50 ) 

C T d - T v 

Note also that T v = 0, whereas f D = 1. 

The universal limiter can be portrayed in the (f c , TA plane 13 . Figure 15 shows the 
constraint boundaries. Note that first-order upwinding (T f — f c ) marginally satisfies the 
limiter constraints. Use of the universal limiter proceeds as follows: 

(i) First compute some (in general, high-order) estimate for the face value, 7} , and 
find the corresponding normalized values of T f and T c . 

(ii) If the point (f c , T f ) satisfies the limiter constraints, proceed to step (iii); if 
not, reset T f to the nearest limiter-constraint at the same T c value. 

(iii) Reconstruct the unnormalized face value 

7} - T, (T„ - T„) * T„ (51) 

(iv) Use this value in combination with second-order diffusion in computing the total 
flux at the CV face; as, for example, in Equation (32) for the left face. 

ULTRA-QUICK Results 

When the universal limiter (for tight resolution and accuracy) is applied to the QUICK 
scheme (giving ULTRA-QUICK), overshoots and undershoots are automatically suppressed 
without additional smearing of the transition region. This is seen in Figures 16 and 17 — 
which should be compared with Figures 11 and 12, and with Figure 2. Note the clean 
monotonic transition in the high-Pe large-a cases, as compared with the unlimited scheme. 
Smooth-region behaviour remains very good, reflecting the uniformly third-order accuracy 
of the basic algorithm. 

Artificial Compression 

Figure 18 portrays a second-order convection scheme in the normalized-variable 
diagram, conforming to universal limiter constraints. The unconstrained portion of the 
scheme (DCB) consists of second-order central-differencing (DC) for f c < 0.5 

? f - \ (1 * T c ) (52) 
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and second-order upwinding (CB) for f > o.5. 



( 53 ) 


The upper constraint boundary (BA) can be interpreted as first-order t/ownwinding 

T f = 1 (54) 

This convection scheme has a tendency to introduce negative artificial diffusion into 
portions of simulated profiles. This can (artificially) enhance resolution of near- 
discontinuities — a phenomenon sometimes known as “artificial compression. ’’ The scheme 
was originally introduced by Roe 15 and named “Ultra-B”; it is related to Roe’s better-known 
“Super-B” scheme, which is also artificially compressive. In terms of simulating near- 
discontinuities, Ultra-B is indeed quite impressive for a second-order scheme. This is clearly 
seen in Figure 19(a) for the infinite-Pe sharp-transition case (a = 100) ; but note the dis- 
tortion of the initially smooth profile (a = 5) in Figure 19(b). As the profile is convected 
downstream, it becomes more and more step-like. This is due to the negative artificial diffu- 
sion inherent in artificially compressive schemes. Similar artificial steepening effects occur 
with the finite-Pe simulations, as well. This is a serious draw-back of artificial-compression 
methods. The phenomenon can be avoided by using higher-order ULTRA-SHARP tech- 
niques, as described in the next section. 

8. ADAPTIVE STENCIL EXPANSION 

Higher-order monotonic resolution of very sharp transitions could be obtained by using 
ULTRA-5th or ULTRA-7th globally. But in most of the flow domain, such high accuracy 
(and concomitant cost) is not called for. It is of interest, from a cost-effectiveness viewpoint 
to construct an algorithm that would use ULTRA-QUICK in “smooth” regions and auto- 
matically branch to a higher-order scheme locally, as needed. The need for the higher-order 
calculation — and correspondingly expanded stencil — can be determined by monitoring 
some suitable “non-smoothness” parameter. One such quantity that comes to mind 
immediately is the local first-difference (proportional to the gradient) across a given CV 
face. For the left face, this would be 

GRADL = T.. - T t _ Uj (55) 

One also needs to detect local changes in gradient; the symmetrically placed average second- 
difference, defined in Equation (37) for the left face, is suitable for this. 

In smooth regions, both GRADL and CRVAVL will lie below certain pre-assigned 
thresholds; in this case, the basic ULTRA-QUICK algorithm is used — this will take care 
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of the bulk of the flow domain since sharp transitions occur in narrow isolated regions, by 
definition. If CRVAVL exceeds the first threshold, THC1 (= 0.1 in the present study), 
the algorithm branches to ULTRA-5th locally, if it also exceeds THC2 (= 0.7), it branches 
further to ULTRA-7th. If GRADL exceeds THG (= 0.35), ULTRA-7th is used immedi- 
ately. Clearly, other threshold strategies could be used; the procedure adopted here has 
evolved through computational experimentation over several test problems. It should be 
noted that the threshold constants are dimensional; i.e., a change in scale, for example, 
would require a corresponding change in threshold values. This problem can be avoided by 
rescaling the threshold constants with respect to an anticipated maximum absolute value of 
the convected variable occurring within the flow-field of interest (in the Smith-Hutton 
problem |T| m « = 2). 

Figures 20 and 21 show results for the ULTRA-3rd/5th/7th scheme described above. 
Clearly, these are highly accurate results, even on the coarsest grid. As seen in the next 
section, the cost is only slightly more than the basic ULTRA-QUICK scheme — but the cost- 
effectiveness (computational efficiency) is greatly enhanced. 

9. OPTIMAL COST-EFFECTIVENESS 

When dealing with higher-order convection schemes, one obvious question that comes to 
mind is: is it better (in terms of total cost) to use a low-order scheme on a very fine grid 
or a higher-order scheme on a coarser grid? Low-order schemes are relatively inexpensive 
per grid-point, but (as seen in the cases shown in this paper) require extremely fine grids for 
reasonable accuracy. On the other hand the added expense (again, per grid-point) of very- 
high-order schemes may not be totally offset by a concomitant coarsening of the grid. To 
be more precise, assume that a desired level of accuracy has been preassigned for a problem 
that has a known exact solution — such as the infinite-Pe Smith-Hutton problem. Take any 
given convection scheme and solve on successively finer and finer grids until the desired 
level of accuracy has been achieved; simultaneously keep note of the CPU time (representing 
cost) at successive grid refinements. Repeat this process with other convection schemes. 
In this way, the cost for a prescribed global accuracy can be assigned to each scheme. 
Alternatively, one could specify an available computational budget and compute the corre- 
sponding accuracy of each scheme as the grid is refined. 

Figure 22 gives the relevant information for the infinite-Pe Smith-Hutton problem with 
a = 100. In part (a) of the figure, the error, given by Equation (31), is plotted versus N , 
on a log-log scale, for first-order upwinding, ULTRA-second-order upwinding (equivalent 
to the Chakravarthy-Osher scheme described by Sweby 16 ), ULTRA-QUICK, ULTRA-5th, 
and ULTRA-7th upwind schemes, together with the ULTRA-3rd/5th/7th upwind scheme. 
Desired accuracy is shown by a dashed line; for each scheme, the corresponding grid refine- 
ment can be found by interpolation. This is cross-plotted onto part (b) of the figure which 
gives CPU-time as a function of N for each scheme. Part (c) of the figure shows the error 
incurred by each method corresponding to a prescribed CPU-time in part (b). Alternatively, 
part (d) of the figure gives the CPU-time for each method corresponding to the prescribed 
error in part (a). 
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From these results, one sees immediately that first-order upwinding is extremely ineffi- 
cient. Among the global higher-order methods, either ULTRA-QUICK or ULTRA-5th is 
seen to be optimal; the 7th-order scheme is significantly less efficient. The most cost- 
effective strategy of all, though, is to use adaptive stencil expansion over the base third-order 
scheme. This is because the wider-stencil (more expensive) higher-order computation is 
automatically used very sparingly — only where needed: in isolated narrow regions 

involving a relatively few number of grid points. For flows involving only relatively smooth 
profiles (such as the a = 5 infinite-Pe case), ULTRA-QUICK is again found to be optimal; 
in this case, the higher-order wider stencil is not called for. 

10. CONCLUSION 

The Smith-Hutton problem is an excellent test of a numerical convection-diffusion scheme, 
especially in the high-convection regime. Strong streamline curvature and rapid local 
variation of the convected variable represent serious challenges to any numerical scheme. 
The availability of exact analytical solutions in the infinite-Pe case is very useful for a 
comparative error analysis. By choosing small values of a, the test-problem can be used 
for simulating smooth-function behaviour, as well. The present formulation of the problem 
uses a staggered grid, interleaving stream function and scalar nodes. Particular attention is 
paid to the outflow boundary condition, assuring ( dT/dy) 0 = 0, consistent with local 
parabolic behaviour. The solution algorithm is based on explicit time-marching until a 
steady state is reached, although ADI tridiagonal solution of the steady equations can also 
be used 13 . 

Exponential-based schemes such as Spalding’s Hybrid 2 , Patankar’s PLDS 3 , or the 
algebraic approximation of Raithby and Schneider 11 , all revert to first-order upwinding for 
convection with modelled physical diffusion (computed but) ignored or suppressed wherever 
the local component grid Peclet number exceeds 2 (Hybrid) or about 6 (for the other 
schemes, including EDS itself). For most flows of practical interest, the grid Peclet number 
is likely to be far greater than 2 or 6 throughout most of the flow-field; under these 
conditions, exponential-based schemes are functioning as first-order upwinding almost every- 
where. The inherent artificial diffusion of such schemes is clearly evident, especially in the 
large-Pe cases. Slow grid-refinement convergence is also observed; this raises serious ques- 
tions regarding grid-refinement claims made in support of exponential-based schemes. Such 
methods should be viewed as of historical interest only — and should not be used for serious 
practical applications. 

Second-, third-, and higher-order upwind methods share a number of similar properties: 
as the order is increased, transition resolution becomes sharper, but overshoots and under- 
shoots become more pronounced, with secondary ripples forming in the case of very-high- 
order schemes. Smooth function and low-Pe performance was seen to be generally very 
good, with error decreasing with order. The two-dimensional third-order (QUICK) scheme 
introduces transverse-curvature terms into the convective fluxes. Other calculations have 
shown that omission of these terms can incur significant error unless the grid is extremely 
fine. As used in this study, higher-order schemes retain the third-order transverse curvature 
terms but omit higher-order transverse and other cross-difference terms; these are costly, 
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algorithmically complex, and seem to have very little effect on the solution. It was also 
found unnecessary to extend diffusion modelling beyond second order. 

The main draw-back of higher-order schemes is the generation of spurious unphysical 
overshoots and undershoots each side of sharp transition regions. This appears to be the 
main reason for a lack of interest in such methods as compared with essentially first-order 
schemes that produce monotonic, albeit extremely artificially diffusive, results. But it is a 
relatively straight-forward task to incorporate monotonizing flux-limiters into higher-order 
schemes, using the concept of a universal limiter, as described in § 7. In terms of locally 
normalized variables, the universal limiter diagram is a simple triangular region with linear 
extensions on each side. When applied to higher-order convective fluxes, the universal 
limiter produces strictly monotonic results without introducing artificial diffusion and 
concomitant numerical spreading of (what should be sharp) transition regions. The tightness 
of the transition resolution increases as the order of the underlying scheme is increased. 
Using a higher-order ULTRA scheme was seen to be a better strategy than relying on arti- 
ficial compression. The negative artificial diffusion inherent in second-order artificial- 
compression methods such as Ultra-B is responsible for extremely tight resolution of near- 
discontinuities; however, as was seen, it tends to distort smooth profiles into a series of 
ramps and plateaus. In a recent paper, Tzanos 17 has also solved the Smith-Hutton problem 
using a third-order convection scheme with a simple limiting strategy essentially equivalent 
to ULTRA-QUICK (but without transverse curvature terms). Tzanos’ paper also gives 
formulas for a variable (adaptive) grid. The results (for a = 10 and Pe = 1000 or 10) are 
very similar to ULTRA-QUICK results for the same parameter values (slight differences are 
due to transverse-curvature terms and different treatment of numerical boundary conditions). 

Among higher-order ULTRA schemes used globally, ULTRA-QUICK and ULTRA-5th 
were seen to be the best schemes in terms of cost-effectiveness: either lowest cost for a 
prescribed accuracy or lowest error for a prescribed cost, as the grid is refined. Adaptive 
stencil expansion — using ULTRA-QUICK as the base scheme and automatically expanding 
the computational stencil to a higher-order ULTRA scheme locally (as needed) — was seen 
to be an extremely cost-effective technique, giving between fifth- and seventh-order accuracy 
for little more cost than that of the underlying third-order scheme. Optimal setting of the 
non-smoothness monitor thresholds requires some experimentation for each new problem; 
but it appears likely that a more general strategy will evolve as experience is gained with this 
new technique. 

The ULTRA-SHARP strategy is ideally suited to steady-state Navier-Stokes calcula- 
tions, as well 18 . If a turbulence model is used, the physics of the model is faithfully 
represented. Very narrow shear-layers can be accurately simulated without fear of artificial 
smearing or oscillation. It is a straight-forward exercise to extend the algorithms to three 
dimensions; and because of the high accuracy obtainable on very coarse grids, reliable three- 
dimensional simulations should soon become practicable for routine engineering calculations. 
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Figure 1 . — Streamline pattern and inlet temperature profile. 
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Figure 3. — Schematic of staggered mesh. Hollow circles are stream- 
function nodes. Solid dots are temperature nodes. Triangles represent 
external pseudo-nodes for temperature. 
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Figure 4.— Detail of typical internal control-volume cell for the temper- 
ature transport equation. 
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Figure 5. — Outflow numerical boundary-condition treat- 
ment. Tq is determined by interpolating a parabola 
through 7 2 and 7, so that {dT/dy) Q = 0. 
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Figure 6. — Comparison of Spalding's (piece-wise linear) Hybrid method, 
Patankar's power-law scheme (triangles), and Raithby-and-Schneider's 
algebraic approximation (circles) with the exponential-differencing scheme 
(light continuous curve). 
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(b) a = 5, % » 0.073. 

Figure 7. — Infinite-Pe results for exponential-based schemes (equivalent to 
first-order upwinding) on a 40x20 mesh. 
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- 1.00 - 0.50 0.00 0.50 1.00 


(c) Pe = 1 0. 

Figure 8. — Inlet and outlet temperature 
profiles for a = 1 00 using PLDS. 


25 










(a) a = 1 00, % = 0.062. 



(b) a = 5,3 = 0.014. 

Figure 9. — Infinite-Pe results for second-order upwinding on a 40x20 mesh. 
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(a) Pe = oo. 



(b) Pe = 500. 



- 1.00 — 0.50 0.00 0.50 1.00 


(c) Pe = 10. 

Figure 12. — Inlet and outlet profiles for 
a = 1 00 using the QUICK scheme. 
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(a) a = 1 00, % — 0.036. 



(b) a = 5, % = 0.003. 

Figure 13. — Infinite-Pe results using fifth-order upwinding on a 40x20 mesh. 
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(a) a = 100, % = 0.034. 



(b) a = 5, % = 0.005. 

Figure 16. — Infinite-Pe results using ULTRA-QUICK on a 40x20 mesh. 
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(a) Pe = ». 



(b) Pe = 500. 



-too - 0.50 0.00 0.50 1.00 


(c) Pe = 10. 

Figure 1 7. — Inlet and outlet profiles for 
a = 100 using the ULTRA-QUICK 
scheme. 
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Figure 18. — The Ultra-B scheme portrayed in the normalized' 
variable diagram. 
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(a) a = 1 00, % = 0.024. 



(b) a = 5, % = 0.004. 


Figure 20. — Infinite-Pe results using the ULTRA-3/5/7 adaptive-stencii- 
expansion method on a 40x20 mesh. 
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CPU-tlme q Error/N 



(a) Error versus grid (c) Error for a prescribed budget, 

refinement. 



N Order 

(b) Cost versus grid (d) Cost for a prescribed accuracy, 

refinement. 

Figure 22. — Computational efficiency of various nonoscillatory schemes for a = 1 00 
and Pe = =°. 
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